Bulletin of Mathematical Biology

Tumors consist of heterogeneous phenotypic cells, such as normoxic cells, which are highly proliferative, and hypoxic cells, which are less proliferative. Their phenotypic switching depends on tumor microenvironmental factors, such as oxygen and nutrient concentrations supplied by local blood vessels. However, during ongoing angiogenesis, the process of sprouting new blood vessels at the tumor site from pre-existing blood vessels, and how this phenotypic switching affects and impacts tumor growth, remains poorly understood. In this article, we formulate a mathematical model to elucidate the crosstalk between vasculature and tumor cellular heterogeneity during tumor progression. The model results show a strong agreement with the experimental data. Our simulation results demonstrate that ongoing angiogenesis increases tumor growth rate. In addition, we observe that the influence of hypoxic cells on phenotypic switching from normoxic to hypoxic is more pronounced than their influence on the transition from hypoxic to normoxic. Furthermore, we perform a global sensitivity analysis using the Sobol's method to assess the importance of the model's parameters. It highlights that the volume at which blood vessels attain half-maximal rate has the maximum effect on the model.

Spatial population structure and seasonality are both central to the spread of many infectious diseases of plants, animals and humans. While seasonal forcing in transmission often plays an important role in epidemiological models of a wide range of infectious disease, and we now have some theoretical understanding of the dynamical impacts of spatial structure, the combined effects of these two ubiquitous processes has not been examined in detail. Here, we develop a novel model to explore the combined influence of spatial structure and temporal variability on disease dynamics. Spatial structure is represented using a lattice-based approach with near-neighbour interactions, while temporal variability is included through regular, seasonal, variation of the transmission rate. We use bifurcation analysis of a pair approximation of the full spatial model to identify the parameter regimes associated with qualitatively distinct dynamical behaviours. The model exhibits a remarkably wide range of complex dynamics, including limit cycles, quasi-periodic cycles, multi-year cycles, chaotic dynamics and bistability between these different states. In particular, complex dynamics occur when reproduction is predominantly local, with the dynamics depending critically on the amplitude of the seasonal transmission rate. We show how high transmission rates, high birth rates and in particular low recovery rates are requirements for complex dynamics. We predict that SI-type disease interactions in plant pathogen systems will show complex dynamics even with relatively global transmission dynamics.

The study of tissue dynamics has been stimulated during the last decades thanks to the use of quantitative descriptions, with the development of several theoretical and computational frameworks, many of them revolving around the notion of reaction-diffusion systems, eventually with additional structure variables beyond time and space. The use of structure variables can accommodate phenotypic traits. In this work, we study a family of competition models, where a given population depends on a resource (e.g. oxygen) and several populations are competing for it. Our quantitative description incorporates phenotypic traits and heterogeneity at the level of cell cycle variations, which influence replication rates via oxygen consumption. This enables us to replicate the fitness of specific subpopulations to environmental conditions (e.g. oxygen shortage or external influences). Using numerical simulations, we show that such models display dynamical pattern formation in the form of coupled travelling wave profiles that expand or retreat at the same wave speed. The full theoretical analysis of such dynamics is quite involved; to circumvent this difficulty, we introduce a quasi-stationary approximation for the resource dynamics. We find that this approximation can reproduce the overall behaviour very accurately, with the additional benefit of allowing theoretical treatment of the reduced model. In this way, we provide estimates on the wave speed which are numerically shown to be robust across a wide range of macroscopic parameters of the full model. The wave speeds are thus found to depend strongly on the proliferation rate of the fittest population, resembling a winner-takes-all dynamics.

After introduction, many non-native marine species are dispersed planktonically. Secondary spread within the non-native range has been shown to prevent the establishment of the introduced species if the advection of larvae prevents sufficient return of larvae to maintain the population in the face of competition with native species. However, those studies have largely neglected the effects of spatial variation in alongshore larval transport. We examine the introduction of a novel species with planktonic dispersal into a more realistic coastal environment which includes spatial variation in larval transport estimated from the Mercator Ocean 1/12th degree global circulation model. The introduction may either be from a distant habitat, or through range expansion. We find that there are locations in the global coastal ocean where introduced species are more likely to persist because of spatial variation of coastal currents. These include regions where alongshore larval transport diverges, such as estuaries. The location where a non-native species is introduced may not be where it flourishes - it cannot be assumed that the region where invading species are first noticed to be abundant is the region where it was introduced. We extend closed-population theory to open coastal systems to estimate persistence as a function of local circulation, habitat extent, and the competitive advantage of the introduced species. Software is provided which allows the estimations of regions where introduced species are more likely to persist and flourish as a function of larval depth behavior, planktonic duration and release timing.

Tumors evolve under constantly changing microenvironmental conditions, and a key major mechanism by which cancer cells adapt to fluctuating environments is stochastic phenotypic switch via epithelial-mesenchymal plasticity (EMP). Here we demonstrate the existence of EMP in human breast epithelial cell lines and assess the fitness of E and M cells in environments that fluctuate between favorable and harsh conditions. We then develop a theoretical framework for analysis of tumor evolution with E-M stochastic phenotypic switch that integrates the behaviors of E and M cells. We characterize the outcomes of tumor evolution via three characteristic times describing i) tumor growth, ii) phenotypic adaptation and iii) fluctuations of the environment. We identify a tradeoff between tumor growth and survival, with different phenotype switching rates maximizing each of these objectives. An anticancer therapeutic strategy is proposed based on tumor survival analysis, revealing that blocking mesenchymal to epithelial transition, rather than epithelial to mesenchymal transition, is critical under fluctuating conditions. Thus, this study elucidates fundamental evolutionary mechanisms of tumor phenotypic adaptation in fluctuating environments with potentially important clinical implications.

Some organisms exhibit declining mortality and increasing fecundity following sexual maturity, a demographic pattern known as negative senescence. According to life history theory, ageing occurs because resources are preferentially allocated to reproduction over somatic maintenance. Models connecting indeterminate growth to negative senescence exist, but none integrate somatic maintenance as a competing allocation decision alongside growth and reproduction. We formulate a life history model in which an individual allocates energy among reproduction, somatic growth, and somatic maintenance and mortality rate depends on both body size and somatic damage. We show that negative actuarial senescence, whereby mortality declines with age, occurs when the proportional change in reproductive value exceeds the proportional change in fitness returns from current investments into reproduction and soma. We derive the necessary conditions for an uninvadable allocation strategy using invasion analysis and Pontryagin's maximum principle, and examine biologically relevant cases numerically. We show that both negative senescence and indeterminate growth arise together as uninvadable outcomes even when maintenance competes for the same resources as growth and reproduction. We show that both diminishing returns to reproduction and diminishing returns to growth can give rise to negative senescence. These results extend the disposable soma theory to organisms with indeterminate growth, in which mortality decreases with size, and identify key mechanisms for the empirically observed association between indeterminate growth and non-senescent demographic trajectories.

Despite the impacts of within-host disease dynamics on disease outcomes in individual hosts and disease spread among-hosts, generic models of within-host population dynamics have received far less attention than their among-host counterparts. While a number of models have been proposed to explore theoretical eco-evolutionary dynamics, they have yet to be evaluated for estimability, raising questions about their ability to provide reliable inference when confronted with data. We evaluated the estimability of two generic within-host population dynamics models by assessing: (1) parameter estimation, our ability to recover correct values of model parameters from data, (2) the consequences of mis-assigning the underlying mechanistic model on parameter estimation, and (3) the reproduction of qualitative dynamics, or, our ability to use parameter estimates to reproduce observed dynamical behaviours. In some cases, fitting a mis-matched mechanistic model to time series data produced reasonable parameter estimates that were able to reproduce system dynamics, and that when provided the data-generating model, parameter uncertainty can produce substantial behavioural uncertainty. Our findings highlight the impacts of structural, parametric, and behavioural uncertainty on inference, and demonstrate the value of improving system-specific knowledge to prevent the use of incorrect functional forms and of measuring consequential parameters to improve estimability.

For organisms with biparental care, successful reproduction hinges on coordination between partners. Seabirds face an extreme coordination challenge because parents must schedule nest attendance on land with long-distance foraging trips at sea. We present a computational model of incubation schedules for a vulnerable seabird, the Leachs Storm-Petrel (Hydrobates leucorhous). Using only simple energetic rules and parameters, the model recapitulates natural incubation rhythms, exposes a tradeoff between parent energy and egg attendance, and predicts severe reproductive failure in harsh environments. Incubation primarily fails through "schedule breakdown" -- a single point in the season when both parents spend too long foraging and the egg dies from cold. The resilience of the developing offspring to neglect is thus a fundamental adaptation to the uncertainties of biparental care. These results raise new alarms about the indirect causes of reproductive failure in sensitive marine species and provide theoretical foundations for the evolutionary ecology of scheduling behaviors.

Antiretroviral therapy (ART) stock-outs interrupt treatment, increase the risk of virologic failure and drug resistance, and erode the population-level benefits of viral suppression. India's National AIDS Control Organization (NACO) manages one of the world's largest public ART programmes, where regimen transitions, evolving formulations, changing treatment guidelines, and procurement-driven fluctuations in drug consumption complicate forecasting. We developed an end-to-end, regimen-specific forecasting workflow to support procurement planning during such periods of instability. We analyzed monthly national ART consumption data from January 2013 through December 2024. A privacy-preserving synthetic dataset was used for pipeline development, followed by final evaluation on real national consumption time series. We compared three model classes, comprising five models: (1) classical models (Holt-Winters and ARIMA), (2) transformer models (TimesFM, which is a large pre-trained time-series foundation model, and its variant with logarithmically transformed values), and (3) hybrid models (variants of a hybrid ARIMA-TimesFM residual model). While the forecast horizon of 18 months remained constant, the train-test period varied across real and synthetic data, as real data was only available until February 2024. For synthetic data, models were trained through June 2023 (test window was July 2023-December 2024), while for real data, models were trained through August 2022 (our test window was September 2022-February 2024). We reported signed percentage deviation to preserve whether models tended to over-or under-predict, and selected models by the smallest absolute deviation. We then derived a regimen-specific model-error buffer, applied only to held-out under-prediction, and deployed the workflow through a no-code dashboard. Forecasting performance was determined using signed percentage deviation (SPD), wherein positive change represents under-prediction and negative change represents over-prediction. Performance varied across regimens, indicating that no single approach was best-performing for all formulations. On synthetic benchmark data, the smallest absolute deviations ranged from 0.46% for adult ABC+3TC to 11.92% for adult AZT+3TC. On real consumption data, classical methods remained competitive for some series, whereas transformer and hybrid models produced better predictive outcomes for others. For instance, for adult AZT+3TC, the Hybrid 70th percentile achieved an SPD of -2.02%, in contrast to the error range of [-15.7, 8.87] for other models. For adult Ritonavir, the ARIMA-TimesFM hybrid at the 30th percentile achieved an SPD of -5.2%, in contrast to the error range of [-14.94, 17.25] for other models. Several formulations, particularly low-volume and transition regimens, nevertheless remained difficult to forecast accurately, underscoring persistent operational uncertainty. This was especially evident across the three pediatric regimens, where all models deviated systematically in the same direction - a more concerning pattern than mere magnitude. For pediatric ABC+3TC, all models over-predicted within a narrow band of [-82.74, -67.43], while for pediatric AZT+3TC and LPV/r 125 mg, all models under-predicted, with ranges of [24.93, 63.73] and [18.32, 52.07] respectively. These findings support a portfolio approach to forecasting in national HIV programmes. Rather than replacing established public-health procurement systems, regimen-specific model selection, directional error reporting, and cautious model-error buffering can strengthen decision support during regimen transitions and other periods of unstable demand.

Neutrophil swarming is a critical immune response in mammals and fish, in which neutrophils are recruited to inflammatory sites where they coordinate into a swarm that neutralizes pathogens. While excessive swarming can drive prolonged inflammation, a quantitative understanding of swarming dynamics remains limited. We developed a one-dimensional radial reaction-diffusion model of neutrophil swarming with two kinetic parameters, in order to capture the self-limiting swarming dynamics in both murine and human neutrophils in response to different inflammatory stimulus sizes. To ensure that the inverse problem is well-posed, we first performed sensitivity and identifiability analyses. We then developed a physics-informed neural network (PINN) to infer the key parameters governing swarm expansion and self-limitation. To account for uncertainty in noisy experimental measurements, we further extended this framework to a Bayesian PINN (B-PINN), which provides credible intervals for the inferred parameters. Both models were validated against synthetic data generated by numerical simulation and subsequently applied to in vitro experimental data from human and murine neutrophils in response to three bioparticle cluster sizes. The PINN-inferred dynamics show that larger bioparticle clusters are associated with greater cumulative recruitment and larger swarms in both species. The models further reveal species-specific differences in both the amplitude of initial recruitment and the timescale on which it self-limits. Additionally, the B-PINN posterior distributions quantify uncertainty in these species- and cluster size-dependent trends and identify where additional measurements would be most informative. To our knowledge, this is the first application of physics-informed machine learning to model neutrophil swarming dynamics. This framework provides a starting point for systematically comparing recruitment dynamics between human and murine neutrophils and offers guidance for future experimental design.

Time compaction theory is a general framework explaining how a brain can efficiently deal with dynamic situations occurring in, e.g., sports games. It involves a geometric representation of the time dimension, which enables effective learning and strategic action planning. The theory has recently received experimental support in humans. However, its current computational model has an important limitation: it does not account for deliberate waiting and speed modulation, behaviors ubiquitous in natural environments. This work substantially extends the original model formulation by a dimensional lifting of an n-D workspace into (n + 1)-D mental space, where time remains geometrically embedded. The proposed biologically inspired computational model can generate adaptive behavior across increasingly complex situations, from navigation in everyday social environments to competitive sports. Furthermore, by actively conditioning the expected responses of other agents and stabilizing future predictions, we introduce the concept of uncertainty points in sequences of generalized cognitive maps to support the generation of adaptive strategies in interactive environments, where future prediction has a limited time horizon. Thus, we provide a mechanism for chaining short-term solutions into long-term strategies, which is illustrated by simulating the behavior of a player in a real football game. Author summaryHumans often anticipate future interactions in dynamic environments. Many behaviors, such as avoiding other pedestrians, letting someone pass through a narrow corridor, or reproducing the kind of dribbling maneuvers performed by elite football players, require deciding not only where to move but also when to move. Existing theories suggest that the brain simplifies such situations by representing future interactions as static spatial maps, making them easier to learn and recall. However, current computational models cannot naturally account for common behaviors such as waiting, slowing down, or modulating speed. Here we show that these behaviors readily emerge if the model space is extended by an additional virtual coordinate that encodes accumulated waiting rather than physical time. The proposed model simultaneously admits a wide variety of behaviors, including speed modulation, multigoal decisions, and compound actions, while preserving the principles of time compaction. We illustrate the model in everyday situations and by reproducing two real football plays, comparing the observed behaviors with model simulations. Our results suggest computational principles through which the human brain may efficiently represent, memorize, and exploit dynamic situations.

Determining the amount of information a neuron receives per unit of time is key to understanding brain connectivity and how neural networks encode and transmit information. In particular, estimating this information flow by distinguishing between excitatory and inhibitory synaptic contributions is critical to understanding neural network function, as maintaining the excitation-inhibition (E/I) balance regulates neuronal excitability and circuit stability, whereas its disruption can lead to a plethora of brain disorders, including neurodegenerative and psychiatric conditions. However, because synaptic conductances cannot be measured directly, inverse methods are required to infer them from the membrane potential --a readily measurable quantity. Although partial solutions have been proposed, accurately estimating these conductances remains a significant challenge due to the complexity and diversity of the inputs. This is particularly true in the spiking regime, where neurons actively fire. In this work, we introduce a novel computational strategy that combines two critical metrics extracted from the time course of the membrane potential recording: the amplitude of the spike and the interspike interval. By using these quantities, the proposed method enables the accurate separation of excitatory and inhibitory contributions, yielding highly favorable results in the spiking regime. Author summaryQuantifying the continuous stream of inputs a neuron receives is key to understanding brain connectivity. Inside the brain, individual cells must maintain a tight balance between excitation and inhibition (E/I) to process information correctly, as any disruption in this equilibrium can impair its functionality. However, directly measuring the underlying excitatory and inhibitory synaptic conductances is technically challenging, and existing mathematical tools often fail when neurons enter their active firing regime. In this work, we introduce a novel computational strategy designed to extract and separate these time-varying conductances directly from the neurons spiking activity. By dynamically tracking just two accessible metrics - the amplitude of the spikes and the time intervals between them - our algorithm estimates both conductance profiles with high precision. Furthermore, we demonstrate that this procedure is highly robust against realistic experimental noise and data variability, providing an accessible framework that does not require complex hardware or an unfeasible number of repetitive experimental trials. By tracking changes in the E/I ratio of the synaptic input, this method provides an efficient approach to detecting pathological imbalances and understanding how local connectivity shapes cellular functionality.

Classical neuronal cable theory relies on quasi-static electric field approximations and neglects magnetic induction, Lorentz force coupling, and transient electromagnetic currents, limiting its ability to fully characterize action potential propagation within geometrically branched axons and dendrites. This work develops a coupled Maxwell-electromagnetic cable framework by integrating finite-difference time-domain (FDTD) solutions of Maxwells equations with extended Hodgkin-Huxley and Fitzhugh-Nagumo membrane dynamics, incorporating magnetic gating perturbations, electromagnetic trans-membrane currents IEM, and nanoscale quantum corrections for thin neural segments. Controlled propagation experiments are designed to quantify deviations from standard cable predictions across asymmetric and symmetric axonal bifurcation geometries. Numerical results demonstrate that inductive magnetic effects lower the critical branch radius for junction conduction failure and break symmetric action potential invasion in geometrically identical child branches under external transverse magnetic fields. An electromagnetic corrected geometric ratio GREM is proposed to revise impedance-matching conditions at branch points, accounting for size-dependent axial current imbalance induced by magnetic and displacement currents. Parent axon conduction velocity deviates substantially from the canonical [Formula] scaling law when electromagnetic feedback and quantum charge distributions are included, triggering early signal blockage at large cable diameters. Collectively, this study establishes that quasi-static cable models underestimate electromagnetic corrections to propagation speed, waveform shape, and bifurcation transmission fidelity; the coupled Maxwell-cable framework provides a comprehensive multi-physics tool for modeling electrodynamic signal behavior in complex neuronal architectures.

A novel mathematical framework to define the threshold of action potentials in excitable cells is presented. Unlike previously applied methods that rely on approximations or bifurcations, the approach focuses on the geometry of membrane potential trajectories. The changes in concavity during the upstroke of an action potential can be directly obtained from a time series of voltages. The concavity criterion is then extended to models based on autonomous dynamical systems where the changes in concavity can be obtained analytically from a curve of inflection points in phase space. The inflection point manifold defines a region required for excitability: all the orbits that cross it contain action potentials, and all the trajectories that contain action potentials are in it. This analytical principle can then be used to define excitability in a dynamical system, and also a measure of excitability that enables quantification and comparisons of excitability across dynamical system. The measure provides a way to compare the excitabilities of systems that model neurons with different electrophysiological phenotypes and consider different stimulus conditions. The traditionally vague physiological concept of electrical excitability is transformed into a rigorous analytical description by considering the time-dependent curvature of the membrane potential. The criterion is robust across smooth, single compartment models of electrical excitability and can be can be extended to single compartment models in higher dimensions, and multicompartment models as well.

High-dimensional population balance equations (PBEs) provide a natural framework for modeling heterogeneous cell populations, but their direct numerical solution becomes computationally prohibitive when the internal state space contains many molecular variables. We propose a hybrid mechanistic-machine learning framework for reducing and simulating PBEs defined over high-dimensional intracellular coordinates. The cell population is described by a number density n(x, t), where x [isin] [R]N represents gene and protein states associated with macrophage activation. A dynamics-preserving autoencoder maps this state space to a low-dimensional latent coordinate z [isin] [R]d, with d << N, while retaining key qualitative features of the underlying gene regulatory network, including attractor structure and multistability. Mechanistic information from the original regulatory dynamics is used to construct interpretable drift and diffusion terms for the reduced latent-space PBE. The reduced PBE is solved using a stochastic Lagrangian particle representation, in which particles evolve according to stochastic differential equations (SDEs) corresponding to the latent drift and diffusion fields. The resulting latent-space solution is subsequently decoded and propagated back into the original state space to recover physically interpretable cellular dynamics. We demonstrate the framework on macrophage polarization under cytokine-dependent regulation, including gene knockout perturbations. Overall, the proposed framework provides a computationally tractable and mechanistically interpretable route for integrating single-cell genomic data with population balance models of cell-state dynamics.

Human Papillomavirus (HPV) infection among men who have sex with men (MSM) has become a significant public health concern, particularly in countries where male vaccination is unavailable. Given the high susceptibility of MSM to HPV and anal cancer, and the unavailability of HPV vaccination for males in low- and middle-income countries (LMICs), there is a need to identify alternative interventions for reducing disease transmission and burden in this population. The novel mathematical model presented in this article couples smoking behavior dynamics with HPV transmission and anal cancer progression among MSM. Smoking reduction is introduced as an intervention to assess its effects on disease transmission and burden. The basic reproduction number (R0) is derived using the next-generation matrix method, and a global sensitivity analysis is performed using partial rank correlation coefficients (PRCC) to identify the influence of model parameters on RR0. Further, the theoretical analysis of the model reveals a backward bifurcation, implying that RR0 < 1 is necessary but not sufficient to eradicate the disease. The study finds that smoking reduction among MSM reduces HPV infection and anal cancer burden relative to baseline projections without intervention. The joint effect of smoking reduction and vaccination shows that the critical vaccination coverage needed to achieve RR0 <1 decreases as the level of smoking reduction increases. A similar outcome is observed for contact reduction. These findings highlight the importance of concurrent interventions, which can significantly curtail the spread of HPV and reduce disease burden in both the high-risk group and the general population.

Tumor spheroid mechanics arise from both the mechanical properties of individual cells and the adhesive interactions that organize them into tissues. The relative contribution of these two factors to the bulk mechanical behavior, however, remains difficult to disentangle experimentally. Here, we develop a computational model of micropipette aspiration to compare the mechanical response of isolated cells and multicellular spheroids within a common computational framework. By independently varying single-cell stiffness and cell-cell adhesion, we quantify their effects on aspiration dynamics, effective elastic modulus, and viscoelastic relaxation. Our results show that increasing single-cell stiffness substantially alters the mechanics of isolated cells but has limited influence on the effective elastic modulus of multicellular spheroids. In contrast, changes in cell-cell adhesion produce pronounced effects on spheroid effective elastic modulus. Nevertheless, both parameters increase the retardation time governing the transition from the initial elastic response to long-time viscous deformation. These findings suggest that multicellular elasticity is governed primarily by intercellular mechanical coupling, whereas the dynamical response to applied stress depends jointly on cell-scale mechanics and cell-cell adhesion.

Insecticide-treated net (ITN) use - defined as the proportion of a population that use ITNs - is a measure of ITN uptake that is used in the estimation of malaria burden and evaluation of intervention programs. However, binary classification of individuals as users or non-users does not account for variations in ITN-related protection attributed to deleterious factors such as chemical and physical degradation, and increased insecticide resistance in vector populations. In this paper, we present a parsimonious model for malaria dynamics in mosquito-human populations in the presence of varying ITN use conditions. Using this model, we propose a new standardised measure of ITN coverage termed the "efficacy-adjusted use" defined as the equivalent level of use, assuming fully efficacious nets, that would be required to achieve the same level of theoretical EIR reduction. This more nuanced measure is used as a proxy for studying ITN-attributed protection across 44 African countries. We find that estimated protection levels in current ITN paradigms is significantly lower than indicated by crude ITN use metrics, with insecticide resistance having the largest deleterious effect. Furthermore, recent adoption of next-generation nets is estimated to have mitigated a 13% reduction in protection compared to a counterfactual pyrethroid only scenario.

A variety of connected systems, ranging from the cytoskeleton to human organizations, dynamically rearrange themselves in order to move through physical or abstract space. However, our understanding of how systems-level behaviors arise from local restructuring actions remains limited, necessitating comparison of real-world data to models that predict network structure and dynamics. To understand these systems, we study an accessible example, the branching slime mold Physarum polycephalum, by imaging the organism as it travels and extracting key fundamental quantities from its continuously remodeling tubular network. By using these quantities as input parameters to a traveling network model, we find that with no further fitting, the model quantitatively matches key emergent properties from P. polycephalum dynamics including path length, relocation time, and search efficiency at different spatial resolutions. These findings demonstrate how a traveling network model can capture P. polycephalum behaviors, highlighting the potential to use traveling networks more broadly for understanding and predicting connected dynamic systems by linking local measurements to emergent, system-wide behaviors.

Auxins are a family of plant hormones involved in various processes across plant tissues and species. The Nuclear Auxin Pathway (NAP) consists of interacting transcription factors (ARFs) and repressors (Aux/IAAs), which govern an individual cells response to changes in auxin concentration. These components are present in all land plants, and many species possess multiple copies of each signalling component. We present a general framework for ODE-based models of NAP submodules with the flexibility to model the promotion and repression of target genes by any combination of transcriptional regulators. We analyse published data and show that auxin treatment in Arabidopsis thaliana roots triggers a range of characteristically distinct temporal response profiles--for both target genes and the signalling components themselves. Using our modelling framework, we recapitulate aspects of this behaviour by presenting examples of real and theoretical NAP subnetworks, and by analysing the effect that these network dynamics have on auxin-mediated transcriptional responses. This work demonstrates the utility of our modelling framework as a general-purpose tool for understanding the function of certain protein-protein and protein-DNA interactions through their effects on the NAP. This exploration of the rich dynamics of more complex signalling pathways promises to advance our understanding of the NAP.